Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 6 December 2008 (MN WF$i style file v2.2) 



Secular evolution of viscous and self-gravitating 
circumstellar discs 

E. I. Vorobyov 1 ' 2 *, Shantanu Basu 3 

OO ' 1 The Institute for Computational Astrophysics, Saint Mary's University, Halifax, NS, B3H 3C3, Canada 
2 Institute of Physics, South Federal University, Rostov-on-Don, 344090, Russia 
, 3 Department of Physics and Astronomy, University of Western Ontario, London, Ontario, N6A 3K7, Canada 

(N ■ 

o : 

. 6 December 2008 

Q 

ABSTRACT 

We add the effect of turbulent viscosity via the a— prescription to models of the 
^ c"| self-consistent formation and evolution of protostellar discs. Our models are non- 

O |. axisymmetric and carried out using the thin-disc approximation. Self-gravity plays 

an important role in the early evolution of a disc, and the later evolution is deter- 
mined by the relative importance of gravitational and viscous torques. In the absence 
of viscous torques, a protostellar disc evolves into a self-regulated state with disk- 
ed | averaged Toomre parameter Q ~ 1.5 — 2.0, non-axisymmetric structure diminishing 

with time, and maximum disc-to-star mass ratio £ = 0.14. We estimate an effective 
viscosity parameter a c g associated with gravitational torques at the inner boundary 
of our simulation to be in the range 10~ 4 — 10~ 3 during the late evolution. Addition of 
viscous torques with a low value a = 10 has little effect on the evolution, structure, 
and accretion properties of the disc, and the self- regulated state is largely preserved. 
\ A sequence of increasing values of a results in the discs becoming more axisymmetric 

in structure, being more gravitationally stable, having greater accretion rates, larger 
sizes, shorter lifetimes, and lower disc-to-star mass ratios. For a = 10 -2 , the model is 
viscous-dominated and the self-regulated state largely disappears by late times. The 
axisymmetry and low surface density of this model may contrast with observations 
! an( i pose problems for planet formation models. The use of a = 0.1 leads to very 

high disc accretion rates and rapid (within 2 Myr) depletion of the disc, and seems 
even less viable observationally. Furthermore, only the non-viscous-dominated mod- 
k>( \ els with low values of a = 10~ 4 — 10~ 3 can account for an early phase of quiescent 

5_j ■ low accretion rate M ~ 10~ 8 Mq yr _1 (interspersed with accretion bursts) that can 

explain the recently observed Very Low luminosity Objects (VeLLOs). We also find 
that a modest increase in disc temperature caused by a stiffer barotropic equation of 
state (7 = 1.67) has little effect on the disc accretion properties averaged over many 
disc orbital periods (~ 10 4 yr), but can substantially influence the instantaneous mass 
accretion rates, particularly in the early embedded phase of disc evolution. 

Key words: accretion, accretion discs - hydrodynamics - instabilities - ISM: clouds 
- stars: formation 
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1 INTRODUCTION 

We have recently demonstrated that disc gravity plays 
an important role not only in the early embedded phase 
of disc evolution but a l so in the late accretion phase 
(|Vorobvov fc Basul 120071 , I2008T ). In the early embedded 
phase, when the infall of matter from the surrounding en- 
velope is substantial, mass is transported inward by the 
gravitational torques from spiral arms that are a manifes- 
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tation of the envelope-induced gravitational instability in 
the disc. In the late accretion phase, when the gas reser- 
voir of the envelope is depleted, the distinct spiral struc- 
ture is replaced by ongoing irregular nonaxisymmetric den- 
sity perturbations. These perturbations produce a residual 
nonzero gravitational torque in the disc. In particular, the 
net gravitational torque in the inner disc tends to be nega- 
tive during first several million years of the evolution, while 
the outer disc has a net positive gravitational torque. This 
is a fundamental property of self-gravitating circumstellar 
discs around low-mass stars. There is also an overall net 
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negative torque in the disc that is related to the removal of 
angular momentum by gas that is accreted in to the central 
object. Although we do not model the gas flow within a cen- 
tral sink of size 5 AU, the angular momentum that is carried 
in to this region may be lost to the system via additional 
processes including magnetic braking and outflows, thereby 
allowing the gas to reach the central star. 

In this paper we seek to determine the effect of other 
mechanisms of radial mass and angular momentum trans- 
port on the secular evolution of self-consistently formed cir- 
cumstellar discs. It is well known that standard collisional 
viscosity (molecular viscosity) is negligible in application to 
circumstellar discs. The best candidate to date is turbu- 
lent vis cosity induced by the m agneto-rotational instability 
(MRI) |Balbus fc H awlcv 1991), though other mechanisms 
such as nonlinear hydrodynamic turbulence cannot be com- 
pletely elim inated due to the la rge Reynolds numbers in- 
volved (e.g. lAfshordi et al.l [20051 ). We make no specific as- 
sumptions about the source of turbulence and parameter- 
ize the magn itude of turbulent viscosit y using the usual a- 
prescription (jShakura fc Sunvaev|[l973l ) 



v — ac 3 H, 



(1) 



where c s is the sound speed and H is the disc scale height. 
The most appropriate value of the parameter a in circum- 
stellar discs is uncertain. Both the shearing box and global 
numerical simulations of the MRI tend to yield the val- 
ues of a that vary significantly from model to model and 
range between 10~ 4 and a few x IP -1 fe.g. lHawlev et ail 
19951; iFleming fc Stonel 120031; iBrandenburg et al.l Il996t 



Stone at al.lll996l : lArmitagelll99Sh . Motivated by these stud- 
ies, we have adopted a similar range of values for a. We also 
assume that a is spatially and temporally constant, i.e. it 
represents a mean value, time-averaged over many orbital 
periods of the disc. Radial variations in a may (and should) 
be present in the disc but it requires a more thorough consid- 
eration of the disc physics and is left for a follow-up paper. 

Our study of the self-consistent formation and long- 
term evolution of circumstellar discs has been p_re- 



ceded by numerical simulations of Lin fc Pringld (ll99cT ) : 
iNakamoto fc Nakagawal l|l995h ; lrIueso fc Guillotl (|2005r i and 
others. However, our work is different in one important as- 
pect - we employ a fully two-dimensional numerical hy- 
drodynamics simulations in the thin-disc approximation (in 
contrast to earlier one-dimensional studies). As a result, 
we account for disc self- gravity self-consistently by solving 
the Poisson integral fsee IVorobvov fc B asu 2006) and need 
not parameterize gravitational torques in terms of the a- 
prescription. Use of the thin-disc approximation is, unfortu- 
nately, a necessity. Fully three-dimensional numerical simu- 
lations are too computationally expensive to study the disc 
evolution on time scales of several Myr. 

The paper is organized as follows. Section [2] gives a 
brief description of model equations and initial conditions. 
The main results are presented in Section [3] The observa- 
tional implications and comparison with earlier studies are 
discussed in Section [5] The main conclusions are summa- 
rized in Section [6] 



2 MODEL DESCRIPTION 

We use the thin-disc approximation to compute the evo- 
lution of non-axisymmetric rotating, gravitationally bound 
cloud cores. We start our numerical integration in the pre- 
stellar phase, which is characterised by a collapsing star- 
less cloud core, and continue into the late accretion phase, 
which is characterised by a protostar/disc system. This en- 
sures a self-consistent formation of circumstellar discs in our 
numerical simulations. Once the disc is formed, its subse- 
quent evolution is determined by an interplay between the 
efficiency of the mass and angular momentum transport in 
the disqj and the infall rate of matter from the surround- 
ing envelope onto the disc. The disc-envelope interaction is 
taken into account self-consistently, since we evolve numer- 
ically the disc and envelope altogether. It means that there 
is no source term in the numerical grid allowing for mass 
deposition from the envelope, but the mass infall rate onto 
the disc is actually determined by the dynamics of gas in the 
envelope. The disc occupies the innermost regions of our nu- 
merical grid and the envelope occupies the rest of the grid. 
We note the infall rate of matter from the envelope onto the 
disc is not necessarily the same as the mass accretion rate 
from the disc onto the protostar. While the former shows a 
fast decline with time, the latter is usually characterized by 
a much slow er decline and has a strong de pendence on the 
stellar mass (|Vorobvov fc Basull2007l . l200Sl) . 

For details of the basic equations, n umerical methods 
and n umerical tests we refer the reader to I Vorobyov fc Basul 
(2006). Here we briefly provide the basic equations modified 
to include the effect of viscosity. The basic equations of mass 
and momentum transport in the thin-disc approximation are 



as 

dt 
,dv p 
~~dt 



= -v P ? + s 9p + (v-n) f 



(2) 
(3) 



where E is the mass surface density, V = J_ z Pdz is the 
vertically integrated form of the gas pressure P, Z is the 
radially and azimuthally varying vertical scale height, v p — 
v r r + v^icj) is the velocity in the disc plane, g p = g r r + 
<?00 is the gravitational acceleration in the disc plane, and 
V p = rd/dr + <pr~ 1 d/d(j) is the gradient along the planar 
coordinates of the disc. We note that g p includes the input 
from the central star when it forms. The viscous stress tensor 
II is expressed as 



n = 2T,u (Vv - ~(V ■ v)e 



(4) 



where Vv is a symmetrized velocity gradient tensor, e is the 
unit tensor, and v is the kinematic viscosity. The compo- 
nents of (V • II)p in polar coordinates (r, (f>) are given in the 
Appendix. We emphasize that we do not take any simplify- 
ing assumptions about the form of the viscous stress tensor, 
apart from those imposed by the adopte d thin-disc approx- 
imation. It can be shown (Lodato 2008) that equation (J3J) 
can be reduced to the usual equation for the conservation of 



1 In fact, discs may also transport angular momentum to the 
external environment due to magnetic braking. This effect will 
be considered in a follow-up paper. 
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angular momentum of a radial annu lus in the axisymmetric 
viscous accretion disc l|Pringlelll98'll ). 

Equations ([2]) and ([3]) are closed with a barotropic equa- 
tion that makes a smooth transition from isothermal to adi- 
abatic evolution at E = E cr = 36.2 g cm -2 : 



fi = 2^o 



V — cfE + c 2 E c 



(e c J 



(5) 



where c s — 0.188 km s _1 is the sound speed in the beginning 
of numerical simulations and 7 = 1.4. Equation ([5]), though 
neglecting detailed cooling and heating processes, was shown 
to reproduce to a first approxi mation the radial tempera ture 
gradients in circumstellar discs l|Vorobvov fc Basul2007l ) and 
the density-temperature relation in collapsing cloud cores 
ijVorobvov fc Basull2006h . 

It should be stressed here that circumstellar discs de- 
scribed by the barotropic equation of state with 7 = 1.4 are 
susceptible to fragmentation and formation of stable clumps 
in the early embedded phase of evolution. A more accurate 
treatment of the energy balance in the disc involving radia- 
tive cooling from the disc surface and shock heating due to 
artificial viscosity has shown that the strength of gravita- 
tional instability in general and the disc propensity to frag- 
mentation in particular depend on the ra t e of cooling (e.g . 



Gammic 2001 



Lodato fc Ricel 



Johnson fc Gammiel 20031 ; iRice et alJl200 j ; 
' 2004 iMei'ia et alJl20oj v In addition, frag- 



mentatio n can be stabilized in the inner discs by slow cool- 
ing (e.g. IStamatellos fc Whitworthl 120081 ) and in the outer 
discs by stellar and e nvelope irradiation (jMatzner fc Levinl 
120051 ; ICai et al.l 2008). However, recent semi-analytical and 
numerical studies (including radiation transfer) reveal clump 
formation in discs (particularly, in their outer parts) around 
stars with mass equal to or more massive than one solar mass 
jKrumholz et al.ll2007l; iMaver et al.ll2007l ; Istamatellos et all 
l2007l : lBosdl200Sl : lKratter et al.ll2008l ). These simulations pro- 
duce discs that are usually hotter than those described by 
a barotropic equation of state with 7 = 1.4 but that does 
not necessarily impl y less efficient tr ansport due to grav- 
itational torques. As ICai et all l)2008h have demonstrated, 
a modest rise in the disc temperature (from the envelope 
irradiation) may in fact promote transport due to a grow- 
ing relative strength of low-order spiral modes (m ^ 2) in 
the disc. In order to examine if a higher disc temperature 
can affect our main results, we consider a stiff er barotropic 
equation of state with 7 = 1.67 in Section [4] 

The kinematic viscosity is computed during numerical 
simulations as 

v = a c s Z, (6) 

where c 2 — dV/dT, is the effective sound speed of (gen- 
erally) non-isothermal gas. The vertical scale height Z is 
determined in each computational cell using an assumption 
of local hydrostatic equilibrium in the gravitational field of 
the central star and the disc (see Appendix A) . 

The initial condi tions are similar to those in 
IVorobvov fc Ba su (2007). The initial radial surface density 
and angular velocity profiles of the model cloud core with 
mass 0.8 Mq and mean molecular weight 2.33 are character- 
istic of a collaps ing axisymmetric magnetically supercritical 
core l|Basulll997l l: 

r S 



(7) Mk) 



- 1 



(8) 



aA 2 + r o 



(7) 



The scale length ro = fcc 2 /(GEo), where k = V2/ir and 
Eo = 0.12 g cm -2 . The central angular velocity is f2o = 
1.1 km s _1 pc " 1 . We have adopted a so mewhat higher value 
of Qo than in IVorobvov fc Basu] ([20071) in order to empha- 
size the burst phase of mass accretion. Our adopted initial 
profiles are characterized by the important dimensionless 
free parameter r\ = f2 2 ro/c 2 . The asymptotic (r 3> ro) ra- 
tio of centrifug al to grav itational acceleration has magni- 
tude V%V (see lBasdfl997h and the centrifugal radius of a 
mass shell initially located at radius r is estimated to be 
r c f = j 2 /(Gin) = v2 V- For our chosen parameters, we find 
77 = 1.42 x 10 , and r f = 16.6 AU for the outermost mass 
shell located initially at r out = 0.04 pc. 

Equations J3), ([5j are solved in polar coordinates 
(r, (f>) on a numerical grid with 128 x 128 points. We use the 
method of finite differences with a time-explicit, operator- 
split solution procedure. Advection is performed using the 
second-order van Leer scheme. The radial points are loga- 
rithmically spaced. The innermost grid point is located at 
r — 5 AU, and the size of the first adjacent cell is 0.3 AU. 
We introduce a "sink cell" at r < 5 AU, which represents 
the central star plus some circumstellar disc material, and 
impose a free inflow inner boundary condition. The outer 
boundary is reflecting. The gravity of a thin disc is computed 
by directly summing the input from each computational cell 
to the total gravitational potential. The convolution theo- 
rem is used to speed up the summation. A small amount of 
artificial viscosity is added to the code, though the associ- 
ated artificial viscosity torques were sho wn to be negligible 
in co mparison with gravitational torques ijVorobvov fc Basul 
l2007fl . A more detailed explanati on of numerical m e thods 
and relevant tests can be found in IVorobvov fc Basul (|2006l . 
l2007l V 



3 RESULTS 

We consider five models, each having identical initial condi- 
tions but distinct values of spatially and temporally uniform 
a. In particular, model 1 is characterized by a = and is 
used as the standard model against which other models with 
non-zero viscosity are compared. Models 2, 3, 4, and 5 have 
a equal to 10 -4 , 10 -3 , 10 -2 , and 10 _1 , respectively. To facil- 
itate the comparison between viscous and non-viscous mod- 
els, a is kept zero in the early evolution in all models and is 
set to its corresponding value only after the disc is formed 
at t » 0.14 Myr. 



3.1 Radial profiles 

Solid lines in Figs [TJ3] show three distinct snapshots in the 
evolution of a circumstellar disc in model 2 (Fig. [1} , model 3 
(Fig. [2]), model 4 (Fig.|3j, and model 5 (Fig. 11} . The numbers 
in the top horizontal row indicate ages of the disc, tdisc- 
Horizontal rows in each figure show (from top to bottom) 
the azimuthally averaged radial profiles of angular velocity 
fi, surface density E, temperature T, Toomre Q-parameter, 
and the relative surface density perturbation AE. Dashed 
lines in each figure give the corresponding radial profiles for 
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Figure 1. Radial structure of the circumstellar disc in the a = 
model 1 (dashed lines) and a = 10 -4 model 2 (solid lines). Top 
to bottom: azimuthally averaged values of SI, £, T, Q, and AS in 
the 0.1-Myr-old disc (left-hand column), 0.85-Myr-old disc (mid- 
dle column) and 1.94-Myr-old disc (right-hand column). The er- 
ror bars in two bottom rows show the minimum and maximum 
A£(r;, <f>i) in each radial annulus. The dotted lines in the sec- 
ond row show the surface density as inferred from the minimum 
mass solar nebula model, and the dot-dashed lines mark a fiducial 
critical density (1.0 g cm -2 ) for transition between the disc and 
envelope. 

the a = model 1. The Toomre parameter is calculated 
as Q = c s f2/(7rG£). In each computational zone (ri,<j)j) we 
calculate the relative azimuthal perturbation to the surface 
density 

M = w^h), (9) 

where N is the number of grid zones in the azimuthal direc- 
tion. The error bars in two bottom rows of Figs lH4l show the 
minimum and maximum AE(ri, (j>i) in each radial annulus. 

It is evident that a = 10~ 4 (Fig. [TJ has little effect 
on the secular evolution of a self-gravitating disc. In the 
end of numerical simulations, when the disc is 1.94-Myr-old, 
the radial profiles of the viscous and non- viscous models are 
nearly identical. The maximum disc radius is about 100 Alfl 
and both the angular velocity and surface density scale as 
r~ 3//2 . In the early evolution (tdisc = 0.1 Myr), the disc is 

2 The disc is defined as the radial distance at which the surface 
density drops below 0.1 g cm -2 



Figure 2. Radial structure of the circumstellar disc in the a = 
model 1 (dashed lines) and a = 10 — 3 model 3 (solid lines). See 
captions to Fig. [l] for details. 

prone to gravitational instability as indicated by both the 
low values of Q = 1.2 — 1.7 and large- amplitude azimuthal 
density perturbations AE = 0.2 — TO. In the late evolution 
(^disc > 0.85 Myr), the disc regulates itself near the bound- 
ary of gravitational stability, Q = 1.7—2.0. This state is 
characterized by ongoing low-amplitude density perturba- 
tions powered by swing amplification at the disc's sharp 
outer edge, resembling truncated circumstellar discs seen in 
young stellar objects. The dotted lines in the surface density 
panels in Figures [T]|4] show the radial surface density profile 
as expected from the minimum mass solar nebula (MMSN) 
mode l, E mmsn = 10 3 (r/AU)~ 3/2 g cm -2 l|Weidenschillind 
Il977h . Our model disc is approximately a factor of ten more 
dense than the MMSN. This implies that most of the disc is 
actually in the optically thick regime, which may have im- 
portant con sequences for the disc mass measurements in T 
Tauri stars (I Andrews fc Williams! 120051 ). 

As the magnitude of a is increased to 10 -3 (Fig. [2}, 
numerical simulations start to show noticeable differences 
between viscous and non-viscous discs in the late evolution 
(idisc <; 0.85 Myr). In particular, the disc starts to spread out 
and its sharp outer edge is replaced with a shallow tail. The 
disc radius at tdisc = 1.94 Myr is twice as large in model 3 
(~ 200 AU) as that in the non- viscous model 1. Both E and 
T decrease throughout the disc. In contrast to non-viscous 
discs, the radial surface density profile in model 3 cannot 
be fitted by a single slope. In particular, E scales as r -0 ' 8 
in the radial range 5 — 60 AU but becomes progressively 
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steeper in the outer portion of the disc. It is also worth 
noting that viscous discs start to slowly drift apart from a 
self-regulation state characterized by a near constant Q ~ 
1.7 — 2.0. For example, the minimum Q value in model 3 at 
tdisc = 1.94 Myr is 2.8. 

As we continue to increase the value of a to 10 -2 
(Fig. [3]), the long-term effect of viscosity on the disc evo- 
lution becomes more profound. It is only the early phases 
of evolution in model 4 (tdi sc = 0.1 Myr) that bear some 
similarities with model 1, the late evolution is considerably 
different. Figure [3] indicates an overall decrease in the val- 
ues of E as compared to those in the non- viscous model 1, 
which results in a disc that is optically thin and cold at 
£disc > 0.85 Myr. The values of E in the radial range 
10 — 200 AU are similar, within a factor of few, to those of 
the MMSN (dotted lines) but become lower in the entire disc 
after 1.0 My r. It is importan t to note that planet formation 
models (e.g. Ilda fc Linlliooj ) require discs with gas surface 
densities a few times greater than that of the MMSN. The 
radial surface density profile cannot be characterized by a 
single slope. It is nearly flat in the inner part of the disc and 
steepens out in the outer parts. The disc radius amounts to 
roughly 350 - 400 AU at t diBC = 1.94 Myr, though we that 
note there is no clear disc boundary and the surface den- 
sity gradually declines with radius to the values typical for 
molecular cloud cores. The disc is virtually axisymmetric in 
the late evolution, except for a small portion near the inner 
boundary, and is characterized by Q 3> 1.0. This marks the 
largest and most noticeable difference between the purely 
self-gravitating disc and the one with a = 0.01 - the latter 
is profoundly gravitationally stable. 

The gravitational stabilization of viscous discs along a 
sequence of increasing a is accompanied by progressively 
more axisymmetric structure. This is an important point to 
keep in mind given the available observational data demon- 
strating that 1-Myr-old discs are non-ax isymmetric and 
show elements of spiral arms and arcs (e.g. iFukagaw a et al.l 
|2004 iGradv et all I2OO1M Furthermore, we point out that 
direct gravitational instability as a scenario for giant planet 
formation is not viable for the viscous-dominated models. 

An increase in a to 10 _1 (Fig. [3} has a catastrophic 
effect on the secular evolution of a circumstellar disc. After 
two million years of evolution, the disc is virtually washed 
out as a result of a very efficient mass and angular momen- 
tum radial transport. In particular, the values of E are more 
than an order of magnitude smaller than those of the MMSN 
at all radii. The Toomre Q-parameter is much larger than 
unity and is off the scale in Fig. [4] at t > 0.85 Myr. Defi- 
nitely, circumstellar discs cannot sustain such large values 
of a for a long time without being completely destroyed. 

Figs I II II show that the angular velocity of the inner 
disc increases with time only by a few per cent. One may 
expect this trend be more profound as the star accumulates 
mass from the disc. However, by the time disc forms in our 
numerical simulations (~ 0.14 Myr), the star has already 
gained most of its final mass (see Fig. [BJ. This is partly 
caused by a relatively low initial rate of rotation of our cloud 



3 Given a large size and relatively low mass of the disc observed 
in AB Aurigae, the nonaxisymmetry in this case may be triggered 
by some kind of close encounter. 
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Figure 3. Radial structure of the circumstellar disc in the a = 
model 1 (dashed lines) and a = 10~ 2 model 4 (solid lines). See 
captions to Fig. [T] for details. 



core and partly by the use of a finite-size (5 AU) sink cell 
in our code. We cannot resolve the very early phases of disc 
evolution. 



3.2 Gravitational torques versus viscous torques 

The differences in the secular evolution of models 1 — 4 
can be understood if we compare the temporal evolution 
of gravitational and viscous torques in these models. The 
blue lines in Figure [5] show the net (positive plus negative) 
gravitational torque in the disc T g , while the red and black 
lines show the net viscous torque T v . The horizontal axis 
shows time elapsed since the disc formation. Both the grav- 
itational and viscous torques are plotted in absolute units 
and the net gravitational torque is always negative. How- 
ever, the net viscous torque is found to change its sign. In 
the early disc evolution at t < 0.1 Myr, the net viscous 
torque is mostly positive and its magnitude is plotted with 
red. In the subsequent evolution, the net viscous torque be- 
comes negative and its magnitude is plotted with black. The 
adopted values of a are indicated in each panel. The net 
gravitational torque T s is found by summing all local grav- 
itational torques r g (r, (j>) — —m(r,(f))d&/d<j) in each com- 
putational cell occupied by the disc, where <& is the local 
gravitational potential and m(r, <j)) is the gas mass in a cell 
with polar coordinates (r, <j>) . The net viscous torque is found 
in a similar manner by summing up all local viscous torques 
T v (r, <f>) — r (V- n)^ S(r, <j>), where S(r, <j>) is the surface area 
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Figure 4. Radial structure of the circumstellar disc in the a = 
model 1 (dashed lines) and a = 10~ 1 model 5 (solid lines). See 
captions to Fig. [l] for details. 



occupied by a cell with polar coordinates (r,<f>). The local 
torques r g (r, <j>) and r v (r, <j>) are actually the local azimuthal 
components of the corresponding forces, acting on a fluid 
element with mass m(r, <j>), multiplied by the radius r. 

Figure[S]clearly demonstrates that gravitational torques 
dominate over viscous torques in the a = 10 -4 model. The 
difference is particularly large in the early evolution and is 
decreasing later but T s remains at least an order greater than 
T v in the 2-Myr-old disc. A qualitative change in the tem- 
poral behaviour of torques is seen in the a = 10 -3 model - 
while the early disc evolution (t < 0.8 Myr) is still controlled 
by gravity, the subsequent evolution is viscosity dominated. 
This effect is even more prominent in the a = 1CP 2 disc, 
in which the viscous torques start to prevail already after 
0.3 Myr of evolution. This explains profound changes in the 
disc structure seen in Fig. [3] The a = 10 _1 model is an 
extreme case, in which viscous torques compete with those 
of gravity even in the early disc evolution and completely 
prevail in the late evolution. 

It is important to note here that gravitational stabil- 
ity properties of astrophysical discs may be modified in 
the pre sence of magnetic fields. For instance, iFromang et al.l 
ll2004al lbl) have demonstrated that the MRI leads to turbu- 
lence, which gives rise to a complicated spiral pattern in 
the disc and lowers the strength of the gravitational stress 
tensor by a factor of 2 due to the nonlinear mode-mode in- 
teraction. In addition, frozen-in magnetic fields tend to de- 
crease the critical Toomre parameter Q CI and the magnetic 
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Figure 5. Temporal evolution of gravitational and viscous net 
torques (in absolute values) in model 2 (upper left), model 3 
(upper right), model 4 (lower left), and model 5 (lower right). 
The horizontal axis shows time elapsed since the disc forma- 
tion. The net gravitational torque (blue lines) is always negative, 
while the net viscous torques can be both positive (red lines) 
and negative (black lines). The torque is measured in units of 
8.66 X 10 40 gm cm 2 s~ 2 . 



disc has to attain a lower v alue of Q in order to destabi- 
lize (|Vorobvov fc Basull2006h . Since a- viscosity is most likely 
caused by the MRI, we acknowledge that the actual strength 
of gravitational torques may be somewhat smaller than that 
shown in Fig.0 We plan to investigate the effect of magnetic 
fields and magnetic braking in a follow-up paper. 

Two important features in Figure [S] need to be empha- 
sized. First, a decline of gravitational torques with time is 
mostly caused by a gradual approach of the disc to a sta- 
ble state. As a consequence, large scale spiral arms in the 
early phase are replaced with on going low-amplitude den sity 
perturbations in the late phase (|Vorobvov fc Basu 200^). As 
Figure[S]demonstrates, turbulent viscosity expedites the on- 
set of gravitational stability and the a = 10~ 2 — 10 _1 discs 
are virtually axisymmetric and gravitationally stable after 
just 1.0 Myr of evolution. Second, a sum of net viscous and 
gravitational torques is always negative. This fact consti- 
tutes a fundamental property of self-gravitating circumstel- 
lar discs, both viscous and non-viscous. The net negative 
torque is related to the ongoing accretion of gas on to a 
protostar - the accreted material removes the disc angular 
momentum, which may be later ejected to the external en- 
vironment via magnetic braking and protostellar jets. As a 
consequence, the rate of change of the disc angular momen- 
tum and the net torque are both negative. It is also inter- 
esting to note that the net viscous torque in the first 10 yr 
of evolution is positive. This is related to the fact that most 
of the disc in this phase is characterized by the dynamic 
viscosity that declines with radius faster than r~ 0,5 . Indeed, 
it can be shown that an axisymmetric Keplerian disc has a 
positive viscous torque, i.e. r (V • II) ± > if fi oc r 13 , where 
< -0.5. 
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3.3 Accretion rates and disc masses 

Figure [5] demonstrates that circumstellar discs evolve 
through two distinct physical regimes during their long-term 
evolution, gravity-dominated and viscosity-dominated. The 
former always precedes the latter and the time when the 
viscosity-dominated regime takes over, or even its existence, 
depends on the value of a. In each of the two regimes, ra- 
dial transport of mass and angular momentum is controlled 
by principally different mechanisms and the mass accretion 
rates on to the central star should bear the imprints of these 
differences. 

The right column in Fig. [6] shows (from top to bot- 
tom) the mass accretion rate in model 1 (a — 0), model 3 
(q = 10~ 3 ), model 4 (a = 1(T 2 ), and model 5 (a = 10" 1 ) 
as a function of time. The mass accretion rate is computed 
as M = — 2%r ivE(r), where v r is the radial velocity at the 
inner inflow boundary r = 5 AU. In all four models the evo- 
lution starts with a sharp rise of M to a maximum value, 
manifesting the formation of a central star at t ~ 0.06 Myr, 
and continues with a short phase of near constant accre- 
tion, in which matter is directly accreted on to the star. 
The disc forms at t ~ 0.14 Myr and the subsequent accre- 
tion history is considerably different between the models. 
The non-viscous model 1 develops short-lived mass accre- 
tion burstfl with M ^ 10 -4 Mq yr -1 , while the quiescent 
phase is characterized by a much lower accretion rate in the 
range a few x 10~ 7 — 10" 8 Mq yr -1 . On the other hand, 
the a = 10 _1 model 5 shows virtually no bursts, while the 
a — 10 -2 model 4 has only a few of them. It is evident that 
the burst activity diminishes along a sequence of increas- 
ing a. There is no quiescent phase of low-rate accretion in 
models 4 and 5 - the mass accretion rate gradually declines 
with time from a few xlO -6 Mq yr" 1 just after the disc 
formation to 10 -8 Mq yr -1 (and lower) at 2.0 Myr. 

There are two observationally valuable consequences of 
the burst phenomenon. First, it can provide an explanation 
for the FU Orionis variables. The apparent lack of bursts in 
the viscous discs will make it more difficult to account for the 
burst activity, though other burst mechan isms may be oper - 
ational in the inner disc at r < 1 AU (e.g. lBeU fc Lm]|l99ll ). 
Second, the existence of the quiescent phase of accretion in 
the early evolutionary phase can potentially be linked with a 
recent detection of Very Low Lumino sity Objects (VeLL Os) 
by the Spitzer Space Telescope (e.g. I Young et al.ll20o3) . A 
previously detecte d object in this category is IRAM 04191 
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Andre et alJll999l ): for another example see lStecklum et al.l 



20071 ). All of these are objects with luminosity L O.ILq 



embedded within dense cores, which, based on previous ob- 
servations, were often classified as starless. Their association 
with dense cores and their low luminosity suggests that they 
are young objects that feature some combination of a sub- 
solar mass and low accretion rate. They might be proto- 
brown dwarfs, but Fig. [6] offers an alternative. The non- 
viscous model 1 can account for VeLLOs as young protostars 
in the quiescent accretion phase with M ~ 10 -8 Mq yr -1 , 



4 The burst phenomenon is caused by disc fragmentation due to 
high rates of mass infall from the enve lope in the early embed- 
ded p hase and is described in detail in lVorobvov fe Basul (|2005| , 
|2006|) . A highly variable ac cretion was also report ed in the case 
of isolated massive discs bv lLodato fc Ricd <2005ft . 
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Figure 6. Right column: temporal evolution of the mass accre- 
tion rate M in the a = model 1 (top), a = model 4 
(middle), and a = 10 _1 model 5 (bottom). Left column: tempo- 
ral evolution of the disc mass (solid lines), stellar mass (dashed 
lines), and envelope mass (dot-dashed lines) in model 1 (top), 
model 4 (middle), and model 5 (bottom). All masses are relative 
to the initial cloud mass M c \ = 0.8 Mq. 



which for the protostellar mass of ~ 1.0 Mq would cor- 
respond to the accretion luminosity 3% that of the solar. 
Clearly, the viscous models 4 and 5 cannot account for VeL- 
LOs, as their mass accretion rates in the early disc evolution 
are too high. The identity of these objects can ultimately be 
revealed by studying t he chemica l evolution of brown dwarfs 
and young protostars (|Leell2007r i. 

Distinct changes in the radial structure of circumstellar 
discs along the sequence of increasing a, as seen in Figs[T} 
21 suggest that viscosity and gravity-dominated discs may 
have quite different masses. The left column in Fig. [6] shows 
disc masses (solid lines), stellar masses (dashed lines) and 
envelope masses (dot-dashed lines), all relative to the initial 
cloud core mass M c \ = 0.8 Mq, in (from top to bottom) 
model 1, model 3, model 4, and model 5. In the a — 
model 1 the relative disc mass reaches a maximum value of 
14% soon after the disc formation and is not significantly 
changing afterwards. On the contrary, disc masses in the 
a = 10 -2 models 4 and a = 10 _1 model 5, though attaining 
similar maximum values soon after the disc formation, show 
a significant decline with time in the subsequent evolution. 
In particular, the disc mass in model 4 drops to only 3% 
that of the central star at 2 Myr, while the disc in model 5 is 
virtually depleted after 2 Myr of evolution. At a first glance, 
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it may seem surprising that viscous and non-viscous models 
have similar (within a factor of 2) disc masses in the early 
evolution although the mass accretion rates in the viscous 
models appear (by eye) to be greater than those in the non- 
viscous model. This paradox is resolved by the fact that the 
bursts are very efficient in regulating the disc mass. During 
each burst a significant amou nt of mass (O.Of — 0.05 Mr) i s 
accreted on to the protostar (| Vorobyov &: Basul l2005. 2006). 

A decline in the disc mass in models 4 and 5 seen in the 
late evolution is caused not only by accretion on to the cen- 
tral star but also by disc expansion due to viscous torques. 
The dot-dashed lines in Fig. [5] illustrate this phenomenon - 
the envelopes in models 4 and 5 are not completely depleted 
by the end of numerical simulations (as in model f and 3) 
but appear to gain some mass in the late evolution. This 
mass is coming from the disc, which dilutes and expands in 
the late evolution (see Figs [3] and [4]| , thus losing part of its 
material to the external environment. The disc expansion 
is explained by the fact that fi in the outer disc scales as 
r 13 , where (3 < —0.5, which causes the disc material to be 
transported outward 

Figure [5] demonstrates an important ingredient of our 
numerical model - our discs are formed self-consistently dur- 
ing numerical simulation and not introduced in the begin- 
ning of numerical simulations. Indeed, the time evolution 
of the mass accretion rate (MAR) in Fig. [6] can be split in 
to three phases. In the first phase, the MAR is character- 
ized by a sharp growth from a negligible value to a peak 
value of approximately 2 x 1O~ 5 M0 yr" 1 . This behaviour is 
characteristic for the runaway collapse phase and stellar core 
formation phase in spherically symmetric cloud core collapse 
simu lations fe.g. lFoster fc Chevalier|[f993l ; I Vorobyov fc Basul 
2005). The second phase is characterized by a near-constant 
MAR at approximately 10~ 5 M© yr _1 . This is typical for 
collapsing self-similar spherically symmetric cloud cores (e.g. 
IShu]|l977| y In these two early phases, the infalling envelope 
accretes directly onto the central star. At approximately 
0.14 Myr a first infalling layer of gas hits the centrifugal 
barrier at r = 5 AU and a centrifugally balanced disc be- 
gins to forrrjf]. Since then, the system enters a qualitatively 
distinct phase of disc accretion. Because the forming disc is 
too small and in the state of near centrifugal balance, the 
MAR drops by many orders of magnitude. This sharp drop 
is a signature of the disc formation and cannot be present in 
numerical simulations if they start right from the disc phase. 
Soon after the disc has formed, it accumulates enough mass 
to trigger radial mass transport and the subsequent time 
behaviour of the MAR is controlled by the interplay of self- 
gravity and viscosity. 



4 STIFFER EQUATION OF STATE 

Numerical simulations with a more realistic treatment of the 
energy balance indicate that the strength of gravitational in- 
stability in general and the disc propensity to fragmentation 
in particular depend on the characteristic ti me of disc cool- 
ing and, hence, on the disc temperature (e.g. lGammie]|200ll ; 



5 In reality, the disc forms earlier but we cannot resolve its evo- 
lution within the inner 5 AU due to the use of the sink cell 
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Figure 7. Radial structure of a purely self-gravitating disc (a = 
0) in the 7 = 1.4 model (dashed lines) and 7 = 1.67 model (solid 
lines). The disc age is indicated in the top frame of each column. 
See captions to Fig.^for details. 



Johnson fe Gammi j 120031 ; iRice et al.ll2003l ; lLodato & R"icl 
2004 iMeiia et al.ll2005l ). In this section we consider the ef- 
fect of increasing the ratio of specific heats from 7 = 1.4 to 
7 = 1.67 in a purely self-gravitating disc. This stiffening of 
the barotropic equation of state should raise disc tempera- 
tures and mimic the effect of a longer characteristic cooling 
time. 

In Fig. [7]we compare the disc radial structure in models 
with 7 = 1.4 (dashed lines) and 7 = 1.67 (solid lines). It is 
evident that an increase in 7 results in a factor of two higher 
gas temperatures in the disc, especially in its inner region. 
At the same time, the gas surface density also increases and 
the resulting Toomre Q-parameter changes insignificantly 
and (in some parts of the disc) may even drop below the 
values characteristic for the colder disc. The radial distri- 
bution of the gas surface density scales as r -1,5 throughout 
most of the evolution, except for the very early phase when 
E oc r~ 2 . The relative azimuthal density perturbations are 
slightly smaller than those in the colder disc. To summa- 
rize, the most substantial change is seen in the radial gas 
temperature distribution. However, since the sound speed is 
proportional to the square root of the gas temperature and 
the gas surface density also increases, the resulting Toomre 
parameter does not change significantly. 

Our next step is to compare the mass accretion rates in 
models with different values of 7. The top and bottom panels 
in Fig. [8] show M in models with 7 = 1.4 and 7 = 1.67, re- 
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Figure 8. Mass accretion rates in the 7 = 1.4 model (top) and 
7 = 1.67 model (bottom) as a function of time. The both models 
are non- viscous (a = 0). The insert shows the rates averaged over 
a period of 2 X 10 4 yr in the 7 = 1.4 model (dashed line) and 
7 = 1.67 model (solid line). 



spectively. A visual inspection of the figure reveals that the 
most noticeable change in the time behaviour of M occurs 
in the early phase of disc evolution between 0.14 Myr and 
0.4 Myr. While the 7 = 1.4 disc exhibits extremely varying 
rates with short bursts followed by longer periods of quies- 
cent accretion, the 7 = 1.67 disc shows a gradually declining 
accretion rate with an order of magnitude fluctuations or 
flickering around mean values. However, the time-averaged 
(over 2 x 10 4 Myr) mass accretion rates are not too dissim- 
ilar, as is seen in the insert to Fig. [7] More specifically, the 
7 = 1.4 disc is characterized by slightly larger /smaller ac- 
cretion in the early/late phase than the 7 = 1.67 disc. It 
is worth mentioning that this tendency is also seen when 
shorter characteristic cooling times are considered. For in- 
stance, the 7 = 1.2 disc[f| has slightly larger/smaller accre- 
tion in the early/late phase than the 7 = 1.4 disc. In any 
case, accretion rates show only a factor of unity difference. 

Why are the time-averaged mass accretion rates rather 
similar in both models? At a first glance, one might expect 
a much large contrast, since an increase in gas tempera- 
ture is expected to moderate gravitational instability and re- 
duce accretion tri ggered by g ravita tional torques. However, 
as was noticed bv lCai et alj (|200&f ). a moderate increase in 
disc temperature may in fact promote accretion due to the 
growing relative strength of lower order spiral modes in the 



6 It is problematic to consider values of 7 smaller than 7 = 1.2, 
since equation J5j tends to overestimate gas pressure when 7 — > 
1.0. 



disc. The global nature of lower order modes makes them 
more efficient agents for the radial mass transport in the 
disc than higher order modes, the latter tend to produce 
more fluctuations and cancellation in the net gravitational 
torque. 

In order to visualize this effect and the strength of spiral 
modes in the disc, we compute the global Fourier amplitudes 
(GFA) defined as 



C n 



£(r, <f>, t) e r dr d<j> 



(10) 



where Ma is the disc mass and rdi sc is the disc's physical 
outer radius. The instantaneous GFA show considerable fluc- 
tuations and we have to time average them over 2 x 10 4 yr in 
order to produce a smooth output. The temporal evolution 
of the time-averaged GFA (log units) for the 7 = 1.4 (top) 
and 7 = 1.67 (bottom) models is shown in the left column 
of Fig. [9] The time behaviour of the GFA is indicative of 
two qualitatively different phases in the disc evolution. In 
the early phase (t < 0.6 Myr), a clear segregation between 
the modes is evident - the lower order mode dominates its 
immediate higher order neighbour in both models. In par- 
ticular, the m — 1 mode is almost always the strongest one. 
The modes also show a clear tendency to decrease in magni- 
tude with time, which explains a gradual decline in the time- 
averaged mass accretion rates shown in the insert of Fig [8] 
The GFA of the 7 = 1.4 disc are somewhat larger than those 
of the 7 = 1.67 disc, except probably for Ci(t), indicating 
that gravitational instability is stronger in the 7 = 1.4 disc. 
The same effect is also seen when we compare the 7 = 1.2 
and 7 = 1.4 discs - the former has larger GFA then the lat- 
ter. In the late phase, however, this clear picture breaks into 
a kaleidoscope of modes, with higher order modes m 2 
competing for dominance with each other and the m = 1 
mode being almost always the weakest one. The GFA have 
dropped in magnitude considerably. However, accretion on 
to the star does not terminate since the effect of these mode 
fluctuati ons is to produce a net negative torque in the in- 
ner disc (|Vorobvov fc Basull2007h . These mode fluctuations 
are not a numerical noise but are rather caused by ongoing 
low-amplitude non-axisymmetric density perturbations sus- 
tained by swing a mplification at the disc's s harp outer edge. 
As was shown by IVorobvov fc Basul l|2007h . self-gravity of 
the disc is essential for these density perturbations to per- 
sist into the late disc evolution. The density perturbations 
(and associated accretion) quickly disappear if self-gravity 
is switched off. 

In the right column of Fig. Owe plot the ratios C* m /Ci as 
a function of time for the 7 = 1.4 disc (top) and 7 = 1.67 disc 
(bottom). It is evident that the relative input from the m ^ 
2 modes is larger in the 7 = 1.4 disc than in the 7 = 1.67 
one. This implies that the mode-to-mode interaction may 
produce more cancellation in the net gravitational torque in 
the colder 7 = 1.4 disc than in the hotter 7 = 1.67 disc. As 
a result, the efficiency of mass transport (due to all modes) 
reduces in the colder disc and becomes similar to that of the 
hotter disc, even though most of the GFA are in fact larger in 
the former than in the latter. This may seem counterintuitive 
but the same effect was fou nd in num e rical s imulations of 
envelope-irradiated discs by ICai et aD (|200Sn . It is worth 
mentioning that a similar modal behaviour was found when 
7 = 1.2 and 7 = 1.4 discs are compared in our numerical 
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Figure 9. Left: Global Fourier amplitudes Cm (in log units) for 
modes m = 1 — 6 as a function of time in the 7 = 1.4 disc (top) 
and 7 = 1.67 disc (bottom). Right: Ratio C m /Ci of the higher 
order modes m Js 2 to the lowest order m = 1 mode for the 
7 = 1.4 disc (top) and 7 = 1.67 disc (bottom). 



simulations - the relative input from the m ^ 2 modes is 
larger in the 7 = 1.2 disc than in the 7 = 1.4 one. The effect 
is not large but is certainly noticeable. 

Alternatively, the similarity of mass accretion rates in 
discs characterized by different 7 may be merely a result 
of the self-regulating nature of embedded accretion discs, 
which re-adjust their gravitational torques (e.g. by increas- 
ing/decreasing surface densities, temperatures, etc.) in order 
to pass on the mass flux coming from the envelope. 

Figure [TO] shows the temporal evolution of viscous and 
gravitational net torques for the same models as in Fig. [5] 
but with 7 = 1.67. The comparison of Figs [5] and [10] indi- 
cates that the time behaviour of the torques in both figures is 
quite similar. There is a slight increase in the strength of the 
viscous torques in models with 7 = 1.67 caused by a higher 
(on averaged) disc temperature than in models with 7 = 1.4. 
In particular, the very early evolution of the a = 0.1 model 
(t ^ 0.02 Myr) is dominated by viscosity, while in the corre- 
sponding model with 7 = 1.4 both the viscous and gravita- 
tional torques are nearly equal. We conclude that a (modest) 
increase in disc temperature does not noticeably affect the 
disc accretion properties averaged over many orbital periods 
(~ 10 4 yr) but can substantially change the instantaneous 
accretion rates. We also note that although the 7 = 1.67 disc 
does not show vigourous bursts of mass accretion (but rather 
an order of magnitude flickering around mean values), the 
bursts can be re-established by increasing the initial ro tation 
velocity of the cloud core (e.g. IVorobvov fc Basull2006l ). 



5 DISCUSSION 

5.1 Constraints on turbulent viscosity 

There are several known physical processes that can con- 
tribute to the radial transport of mass and angular mo- 
mentum in circumstellar discs. These transport mechanisms 
include gravitational torques, either due to internal (self- 
gravity) or external (companion star) sources, and turbu- 
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Figure 10. The same as Fig. [5]only for models with 7 = 1.67. 



lent viscosity. The effect of self- gravity on the secular evo- 
lution of a circ umstellar disc has been considered in our 
previous paper l|Vorobvov fc Basull2007h . In this paper we 
added the effect of turbulent viscosity due to yet unspeci- 
fied source of turbulence. Physical processes in circumstellar 
discs that can, under specific conditions, give rise to turbu- 
lence include the magneto-rotational instability and vertical 
convection. The latter may be heavily suppre ssed in the ex- 
ternally heated discs l|Ruden fc Pollacklll99lJ) . which leaves 
the MRI as the most probable source of turbulence in cir- 
cumstellar discs. lAfshordi et al.l (|2005l ) have recently argued 
on analytical grounds that hydrodynamic turbulence (not 
related to the MRI) can arise in cold Keplerian discs char- 
acterized by Reynolds numbers > 10 4 and having both the 
fine-tuned conditions and appropriate feedback mechanism, 
though this idea has not been proven so far by numerical 
simulations. 

Many numerical efforts to characterize the MRI- 
induced turbulent viscosity qualitatively have been made 
the recent years. These include l ocal shearing box 



models (e.j 



Brandenburg et all Il996l; Istone at al.lll996l ; 
iFleming fc Stond 120031 : ISato et alJ 120041 ) and global mag- 
netohydrodyna mic simulations of accretion discs (e.g. 
lArmitage|[l998l ). though the latter are usually limited by 
a small number of orbital periods. Most authors calculate 
the mean ratios of the Maxwell and Reynolds stresses ver- 
sus midplane gas pressure, T^/Po and T^f/Po, respectively, 
which, to a factor of unity, are proportional to a in Keplerian 
discfl 

The derived valu es of a vary significantl y between the 
studies. For instance, lHartmann et al.l (| 19981 ) derived a « 
0.01 using observed disc sizes and a simple model for the evo- 
lution of an axisymmetric visc ous disc, in which v iscosi ty is a 
power-law function of radius. IFleming fc Stond (|2003l ) have 
employed a stratified model of accretion discs in a shear- 
ing box approximation, in which the upper layers are MRI- 
active while the central regions are quiescent. They found 



7 In fact, this is true only for axisymmetric discs, in which 
= (rilog f2/eilog r)aPo and the rate of strain d log H/ri log r is 
constant (e.g. -3/2 in Keplerian discs). In non-axisymmetric discs 
the relation between and Pq is more complex. 
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that T^ y /Po has a mean value of about 10~ 3 in the MRI- 
active upper layers but drops to a negligible value in the 
midplane. On the other hand, T^f/Po in the vertical direc- 
tion is roughly constant at a few xlO -4 . This implies mean 
values of a in the range 1CT 4 - 10~ 3 . Low mean values of a 
were also reported by Brandenburg et al. (a = 0.007) and 
Stone et al. (a < 0.01). On the other hand, global numerical 
simulations tend to yield larger values for a. For instance, 
lArmitagel (| 19981 ) found mean values between 0.05 and 0.1. 

Large variations in a, both in space and time, imply 
that the development of the magneto-rotational instability is 
strongly dependent on the local conditions in the disc. Nev- 
ertheless, we can still learn about circumstellar discs from a 
simple model used in our paper, if we assume that the con- 
stant a represents a mean value, time-averaged over many 
orbital periods of the disc. Radial variations in a may (and 
should) be present in the disc but it requires a more thor- 
ough consideration of the disc physics (i.e. the ionization 
balance) and is left for a follow-up paper. 

Our numerical simulations unambiguously demonstrate 
that circumstellar discs cannot sustain turbulent viscosity 
with a spatially and temporally averaged a > 0.1. Such discs 
would have vanished during just one million year of evolu- 
tion. The ubiquitous presence of older discs makes such large 
values of a unlikely. On the other hand, low values of a of 
order 10 -4 make little effect on the secular disc evolution, 
which in this case is completely governed by gravitational 
torques rather than viscous ones. In the case of a = 10 -3 , 
viscosity does have some effect on the disc radial structure 
but the magnitude of these changes are modest - the surface 
density profile becomes somewhat shallower in the inner disc 
and at the disc's outer boundary, and the disc size increases 
by a factor of 2 as compared to that of the non-viscous one. 
The a = 0.01 disc sees considerable changes in its radial 
structure in the late evolution, and the gravitational stabi- 
lization of such discs presents difficulty to account for non- 
axisymmetric structure and poses problems for the theoret- 
ical idea of giant planet formation via direct gravitational 
instability. Moreover, the gas surface density in the entire 
a — 0.01 disc becomes lower than that of the MMSN af- 
ter 1.0 Myr. This also poses problems for planet formation 
models, which often require discs with gas surfac e densities 
a few times greater than that of the MMSN (e.g. Ilda fc Linl 
120041 ). We conclude that the mean value of a (averaged over 
many orbital periods) should lie in the range 10" 3 - 10~ 2 , 
although large transient variations around these values can 
still be present in real discs. 



5.2 Effective viscosity due to gravitational torques 

An interesting way to account for the combined effect of 
gravitational and viscous torques is to calculate an 'effective 
alpha' a c ff near the inner boundary of our simulatiorjf]. If we 
apply the steady-state mass accretion rate formula for thin 
viscous discs and also apply the a— prescription, then 



M = 3vn/E = 3na cB c a ZY,. 



(11) 



We calculate M in each simulation at the boundary of the 
sink cell (r = 5 AU), and c s ,Z, and E at some distance 

8 For various methods to define ct e ff see e.g. lLodatol J2007T> 
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Figure 11. Temporal evolution of the effective viscosity parame- 
ter a e ff that accounts for both viscous and gravitational torques, 
for each of our models. 



(r ~ 9 AU) since E decreases somewhat near the sink cell. 
These numbers are used to generate a c a, which is plotted 
versus time in Fig. [TT] for each of our models. The values 
of a e ff may be scaled down by a factor ~ 2 if use also the 
values of M at 9 AU instead of 5 AU. Clearly, gravitational 
torques in the absence of turbulent viscosity accounts for an 
a e a in the range 10~ 4 — 10~ 3 during the late evolution. This 
is why the addition of a viscous a of at least 10 -3 is required 
to see significant changes in the disc evolution. 

It is also interesting to compare our results with previ- 
ous numerical simulations of the secular evolution of v iscous 
circumstellar discs. For instance, iLin fc Pringkl (jl990h have 
considered the formation and evolution of a circumstellar 
disc formed during the collapse of a rotating cloud core with 
initial mass 1.0 M®. They use the usual diffusion equation 
describing the evolution of t he surface density in a viscou s 
axisymmetric accretion disc l|Lvnden-Bell fc Pringleiri974l ) 



9E 
dt 



ld_ 

r dr 



d 



{r 2 Q)' dr 



VEr 3 fi') 



(12) 



complemented by some form of the energy equation describ- 
ing the internal energy balance in the disc due to viscous 
heating, energy input from the accretion process, and disc 
radiative cooling. Our model, though neglecting detailed 
thermodynamics, accounts for a possible disc asymmetry by 
directly solving the corresponding fluid dynamics equations 
for a thin disc. The effective kinematic viscosity in Lin & 
Pringle's model comes from turbulent viscosity and gravita- 
tional instability, the former is parameterized using a usual 
Shakura & Sunyaev a-prescrip tion (eq.[T]), while th e latter 
is taken into account following ILin fc Pringlei fl987): 



2/i 
3 





(Ql 



(#) if Q^Qc 

otherwise. 



(13) 



We need to parameterize only the viscous torques, the 
effect of gravitational torques is taken into account self- 
consistently. 

Both approaches yield circumstellar discs that share 
some common characteristics. For instance, Lin & Pringle's 
model Al (a, n = 0.01) produces a rather cold disc (T ~ 
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10 K) that features a near-flat surface density profile near 
the inner boundary and scales roughly as r -1 at 100 — 
1000 AU. Their disc has a sharp outer edge upon formation 
but it spreads out in the course of evolution. These features 
are also seen in our modeling, though our disc in the cor- 
responding model 3 has a smaller size (~ 300 AU) and a 
lower surface density. This brings about the most striking 
difference fo und between the two approaches - the derived 
disc masses. iLin fc Pringle! (ll990T ) have reported disc masses 
that are an order of magnitude larger than the correspond- 
ing stellar masses in the early evolution. Although this large 
difference reduces with time, the disc mass is still compara- 
ble to that of the star in the late evolution (t > 1.0 Myr), 
irrespective of the value of a. Such massive discs are not 
observed. Our models predict maximum disc-to-star mass 
ratios of £ = 14% (model 1), and this value quickly reduces 
with time for the viscosity-dominated discs (e.g. model 3). 
Although our obtained values of £ for non-viscous models 
still seem to be greater t han those usually inferred for T 
Tauri stars, 0.5%-l% fe.g. lAndrews fc Williamsll2005l ), they 
are not unfeasible given that the measured disc masses may 
be underestimated by conv entional methods by a s much as 
an order of magnitude (e.g. lHartmann et alj|2006h . Figure [5] 
indicates that gravitational torques are a dominant mecha- 
nism of mass and angular momentum transport in the disc 

in its early evoluti onary phase. 

The fact that ILin fc Pringle! <|l990h obtain overmassive 
discs in their numerical simulations suggests that either the 
parameterization given by equation (|13[) is inadequate, cer- 
tainly in the early disc evolution, or the ratio of rotational- 
to-gravitational energies f3 in their model is too large, re- 
sulting in most of the initial cloud mass landing on to the 
disc and through the disc on to the star rather than directly 
on to the star. In other words, the phase of near constant 
accretion of matter from the envelope directly on to the star 
(see Fig. |SJ) is very short in the Lin & Pringle model and 
the disc is not capable of processing the infalling mass to 
the star fast enough to keep its mass low. Indeed, Lin & 
Pringle have adopted /3 in the range 0.25 — 0.64, which is 
much larger t han the values recen tly inferred for molecular 
cloud cores bv lCaselli etaH (|2002h . p = 10" 4 - 0.07. In our 
model, P is set to 1.4 x 10~ 3 . 



6 SUMMARY 

Using numerical hydrodynamics simulations we have stud- 
ied the long-term evolution (2.0 Myr) of self-consistently 
formed, self-gravitating circumstellar discs that are subject 
to turbulent viscosity. We seek to determine the effect of 
viscosity on the radial structure and accretion properties of 
self-gravitating discs around low-mass (~ 0.7 Mq) proto- 
stars. We make no specific assumptions about the source of 
turbulence in circumstellar discs and parameterize the mag- 
nitude of turbulent viscosit y using the usual a-prescription 
|Shakura fc Sunvaevl Il973h . Four models with a spatially 
and temporally uniform values of a — 10~ 4 , 10 -3 , 10~ 2 , 
and 10 _1 were considered and compared with the standard 
model characterized by a = 0. We find the following. 

(i) Low values of a of order 10 -4 make little effect on 
the secular evolution of a self-gravitating disc, the radial 
structure and accretion properties of which in this case are 



completely determined by gravitational torques rather than 
by viscous ones. 

(ii) At values a of order 10 -2 , the discs see considerable 
changes in their radial structure, with a surface density that 
is axusymmetric and has values that are already below that 
of the MMSN by 1 Myr. This is problematical for planet 
formation. 

(iii) High values of a of order 10 _1 make a catastrophic 
effect on the disc secular evolution. Most of the disc mass is 
quickly accreted on to the protostar and the rest is dispersed 
to the external environment. The disc mean surface density 
drops below 1.0 g cm -2 during just 1 — 2 Myr. 

(iv) The net viscous torque in the disc is positive in the 
early evolution (< 0.1 Myr) and negative afterwards. On the 
other hand, the total (viscous plus gravitational) net torque 
is always negative, which is related to the removal of angular 
momentum from the disc by gas that is accreted in to the 
central region. 

(v) Use of a stiffer barotropic equation of state (7 = 1.67 
instead of 7 = 1.4) and associated rise in disc temperature 
can substantially affect the instantaneous accretion rates 
(particularly in the early evolution) but have little effect 
on the disc accretion properties averaged over many disc 
orbital periods (~ 10 4 yr). This is because a decrease in 
the intensity of mass accretion bursts in the hotter disc is 
compensated by an increase in the relative strength of lower 
order spiral modes (m ^ 2), which are the most efficient 
agents for radial mass transport in the disc. 

(vi) Viscous-dominated models have difficulty to account 
for some physical properties of circumstellar discs. For in- 
stance, they become virtually axisymmetric and gravita- 
tionally stable after just 1.0 Myr of evolution. Moreover, 
the lack of a quiescent phase of low-rate mass accretion 
(M ~ 10" 8 Mq yr 1 ) in the early evolution of viscous 
discs will make it difficult to account for Very Low Lu- 
minosity Objects (VeLLOs), which are presumably young 
objects that feature s ome combination of a sub-solar mass 
and low accretion rat e (|Young et al ] |2004l ; lAndre" et al 1ll999l ; 
IStecklum et alj|2007t ). On the other hand, non-viscous self- 
gravitating models can naturally account for both the ap- 
parent disc non-axisymmetry and a phase of very low lumi- 
nosity. 

We emphasize that the a-parameter in our models does 
not include a (possible) contribution from Reynolds stresses 
due to self-gravity. This fact should be taken into account 
when comparing our predicted values of a with those de- 
rived in other studies, which may include a contribution from 
the gravitational Reynolds stresses (e.g. lGammiell200lh . We 
have also not considered a possible contribution of Reynolds 
stresses in th e transport of mass a nd angular momentum. 
According to lLodato fc Rice] (2004), this contribution may 
be comparable to that of the gravitational torques, which 
means that the importance of self-gravity may be underes- 
timated by a factor of 2. 
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APPENDIX A: DISC SCALE HEIGHT 

We derive the disc vertical scale height Z at each computa- 
tional cell via the equation of local vertical pressure balance 



pc s = 2 



z ,gas + ffz.st) dz, 



(Al) 



where p is the gas volume density, gr z , ga s and g ZzSt are the ver- 
tical gravitational accelerations due to self-gravity of a gas 
layer and gravitational pull of a central star, respectively. 
Assuming that p is a slowly varying function of vertical dis- 
tance 2 between z = (midplane) and z = Z (i.e. E = 2 Z p) 
and using the Gauss theorem, one can show that 



P9z, gas dz 



pffz.st dz = 



GM„p 



1 - 



1 + 



S 

2pr 



-1/2 



(A2) 



■(A3) 



where r is the radial distance and M t is the mass of the 
central star. Substituting equations (|A2|I and (|A3|) back into 
equation (|A1|) we obtain 



pet = — GE + 



2GM*p 



1 + 



E 

2p~? 



-1/2' 



(A4) 



This can be solved for p given the model's known c^, E, 
and M* , using Newton- Raphson iteration. The vertical scale 
height is finally derived as Z — E/(2p). 



APPENDIX B: DIVERGENCE OF THE 
VISCOUS STRESS TENSOR 



The components of V ■ n in polar coordinates (r, (j>) are 

(Bl) 



(V-H) r = l± r ILrr + —**<■+- — 
r or r o<p r 



(B2) 



where we have neglected the contribution from off-diagonal 
components H r z and n^ z . 
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